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ABSTRACT 

We obtain equilibrium solutions for rotating compact stars including the special rela- 
tivistic effects. The gravity is assumed to be Newtonian, but we used the active mass 
density, which takes into account all the energies such as motions of the fluids, in- 
ternal energy, pressure energy in addition to the rest mass energy, in computing the 
gravitational potential using Poisson's equation. Such a treatment could be applica- 
ble to the neutron stars with relativistic motions or relativistic equation of state. We 
applied the Hachisu's self-consistent field (SCF) method to find spheroidal as well 
as toroidal sequences of equilibrium solutions. Our solutions show better agreement 
than Newtonian relativistic hydrodynamic approach that does not take into account 
the active mass, with general relativistic solutions. The physical quantities such as 
the peak density, equatorial radii of our solutions agree with general relativistic ones 
within 5%. Therefore our approach can be a simple alternative to the fully relativistic 
one when large number of model calculations are necessary as it requires much less 
computational resources. 

Key words: gravitation, hydrodynamics, stars: neutron 



1 INTRODUCTION 

Finding an equilibrium solution is a starting point for the 
studies of the dynamical evolution of any objects. The se- 
quence of equilibrium models as a function of particular pa- 
rameter could also give us some guess and insight for the 
real evolutionary dynamics when direct numerical simula- 
tions are difficult. It is rather trivial to get solutions for 
non-rotating spherical stars because it only requires the inte- 
gration of a single ordinary differential equation with proper 
boundary conditions if the equation of state is barotropic. 
For rotating stars having non-spherical shapes, however, 
solving the elliptic equation for the problem is not so simple 
because the location of boundary of the star is not prede- 
termined on which the boundary condition should be im- 
posed. Since James (1964) first found the solution of rotating 
stars by directly integrating the elliptic equation, many peo- 
ple have tried to develop more efficient and accurate ways. 
Self-consistent field (SCF) method (Ostriker & Mark, 1968), 
which is one of the possible approaches, uses integral repre- 
sentation rather than solving the differential equations di- 
rectly. Based on SCF approach, Hachisu (1986a,b) developed 
a very successful method which covers almost all possible 
configurations of rotating stars with barotropic equation of 
state and Newtonian gravity. Soon after it was devised, this 
method was adopted to generate general relativistic initial 



data of rotating stars (Komatsu et al. 1989a,b, KEH here- 
after). It was also used for the study of supramassive stars 
which cannot be treated with Newtonian gravity (Cook et 
al., 1992, 1994a,b) and for the study of rotating stars with 
realistic equation of state (Cook et al., 1994a, b; Stergioulas 
& Friedman, 1995). Now, SCF is one of the most popu- 
lar methods to find equilibrium solutions for wide ranges of 
problems other than single rotating stars e.g., binary stars 
(Baumgarte et al., 1998), magnetized stars (Kiuchi & Ko- 
take, 2008). 

Relativistic motions of self-gravitating objects are pos- 
sible in various astrophysical circumstances. During the core 
collapse of proto-neutron stars, the fluid velocity can reach 
~ 0.2c where c is the speed of light. Recently found millisec- 
ond pulsar XTE J1789-285 has rotation frequency of 1122 
Hz (Kaaret et al., 2007), which corresponds to the rotation 
speed ~ 0.2 — 0.4c at the stellar surface. Many phenomeno- 
logical models for neutron stars require relativistic equation 
of state whose pressure and internal energy density can be 
comparable to the rest mass energy density. All these in- 
gredients should be well implemented in general relativistic 
hydrodynamics. It is now possible to study these objects 
in three-dimensional numerical simulations by taking into 
account the dynamics of spacetime [see for a review, Font 
(2008)]. However, it is very difficult to simulate realistic de- 
tails of neutron star physics general relativistically because it 
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requires too much computational resource. Therefore, New- 
tonian gravity is still widely used in cases the gravitational 
fields are not so strong (Tohline & Hachisu, 1990; Ou & 
Tohline, 2006; Petroff & Horatschek, 2008). Sometimes New- 
tonian approach has been adopted even for the problems in 
which relativistic motions or relativistic equation of states 
are involved because of its simplicity. This approach might 
be justified if it does not produce significantly different re- 
sults from the fully relativistic calculations. However, it is 
better to find a new approach which can account for the 
relativistic properties in a generic way. 

The purpose of this paper is to implement the special 
relativistic effects properly in finding the equilibrium solu- 
tions for rapidly rotating compact stars under the Newto- 
nian gravity. Our approach is different from previous stud- 
ies in the sense that we use the active mass which takes 
into account the contributions from various energies to the 
gravitational potential, in addition to the rest mass energy. 
Our framework is based on the Hachisu's SCF method. The 
implememtation of active mass should provide an improved 
agreement with the general relativistic solutions since the 
concept of active mass essentially comes from the general rel- 
ativity. The rotating neutron stars and strange stars would 
be the best examples to which our method can be applied. 

This paper is organized as follows. In section 2, we de- 
scribe our approach taken to find a hydrostatic equilibrium 
solution of rotating stars. The details of numerical scheme 
to find the solutions are given in section 3. We compare our 
solutions with those obtained by Newtonian and general rel- 
ativistic approaches in section 4. The properties of our so- 
lutions for the rotating stars are presented and discussed 
in §5. The final section provides summary and concluding 
remarks. 



2 FORMULATION 

The general relativity is the only way to incorporate the spe- 
cial relativity with gravity. The dynamics of the gravity (or 
spacetime) can be analyzed by solving the Einstein equa- 
tions. At the same time, the equations of motion for matter 
are given by the conservation of the energy-momentum ten- 
sor which is the source of the gravity. All these are highly 
coupled and nonlinear equations so that the state of the arts 
numerical techniques and massive computational resource 
are required. If the gravity is weak, however, one can take the 
so called "weak field approximation" by linearizing the Ein- 
stein equations and solving only the matter dynamics in flat 
spacetime. Although it is now much simpler and sometimes 
is possible to treat analytically, we still have 6 unknowns 
for the spacetime (after gauge fixing) in general and have 
4 unknowns for axisymmetric rotating stars of our concern. 
But this is too many for our purpose. Therefore, just like the 
Newtonian approach, we take only one dynamical variable 
(gravitational potential) for the spacetime by assuming the 
following metric, 



ds 2 = -(1 + 2$)di 2 + (1 + 2<S>y 1 8 lJ dx l dx J , 



(1) 



where $ is the Newtonian gravitational potential. In 
the weak field approximation, the assumption results in 
Tij, Tio <C Too- This condition can be violated by the rela- 
tivistic problems of our concern. Since we try to develop an 



effective method, however, we shall not stick to this discrep- 
ancy. Instead, we take into account the possible relativistic 
circumstances of Tij, To ~ Too by introducing "active mass 
density" , p ac twc which was first suggested by Tolman (1934). 
The opposite concept, "passive mass density is just the in- 
ertial rest mass density which appears as a source term in 
the Poisson's equation. However, it is clear from the general 
relativity and the mass-energy equivalence that all kinds of 
energies can contribute to gravitational potential that influ- 
ences spacetime geometry. The corresponding active mass 
density given by Tolman (1934) has the form 



TryrjiO rr~ii rrtO 

, -^-t () — J-i J-Q 



(2) 



where T is the trace of the energy-momentum tensor. Ac- 
cordingly, we modify the Poisson's equation to have the 
Pactivc in the source term, 



V 2 $ = 47rGp acti 



(3) 



To describe the rotating star we assume that the energy- 
momentum tensor of matter is that of perfect fluid as given 
by 



(4) 



where po is the rest mass density, which is proportional to 
the number density of baryon of the fluid, P is the pressure, 
tt M is the four velocity of a fluid element with respect to the 
Eulerian observer and </ M „ is the spacetime metric assumed 
in eq. (1), and H is the specific enthalpy which is defined as 



H = l + e + 



Po 



(5) 



with e being the specific internal energy. The equation of 
state is assumed to be barotropic, i.e., P = P(po). Further- 
more, for convenience, we only consider simple Polytropic 
equation of state, 



P = Kp Q 



(6) 



where k and N are polytropic coefficient and index, respec- 
tively. 

From the general relativistic version of hydrostatic 
equation found by KEH, we can easily get the hydrostatic 
equation for our metric as follows, 



PoH 



VP + V In (1 + 2$) 



1 - v 2 



Vv + 



v 2 vfi 
i - v 2 IT 



o, 



(7) 



where fi is the rotational angular speed, and v = 
fir sin 6>/ (1 + 2$). The fluid four- velocity u M due to the ro- 
tation is thus given by 



^(1 + 2$) (1 - v 2 ) 



(l,0,0,fi). 



(8) 



The axial symmetry is assumed and the spherical polar co- 
ordinates (r, 9, <j)) are used. If the last term of the left hand 
side of equation (7) is a function of fi, the hydrostatic equa- 
tion is integrable. Usually, the fi-dependence is assumed to 
have the following form (KEH), 



F (fi) = A 2 (fi c - fi) 



fir 2 sin 2 9 



(l + 2$) 2 -fi 2 r 2 sin 2 6> 



(9) 
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where Qc is the angular speed at the rotation axis, and A is a 
constant which controls the behavior of Q. This choice of ro- 
tation law automatically holds the rough Rayleigh stability 
condition i.e., the specific angular momentum (j) should not 
become smaller outwards. Two limiting cases are ^4^00 
and A — > 0. The former choice gives the uniform rotation 
and the latter one leads to rotations with a constant angular 
momentum (KEH). 

By integrating equation (7) with the rotation law given 
in equation (9), we get the integral representation of the 
hydrostatic equation as follows, 



InH + i In (1 + 2$) 



+ Iin(i-« 2 ) -\a 2 (q-q c ) 2 = c, 



(10) 



where C is an integration constant. By imposing boundary 
conditions, this equation can be used to obtain the rest mass 
density. From equations (1), (2), (4) and (8) the active mass 
density is given by 



Pactivc = PqH 



1-v 2 



+ IP. 



(11) 



The gravitational potential is found by solving the modified 
Poisson's equation [eq. (3)]. 



3 METHOD 

3.1 Hachisu's Self-consistent Field 

In addition to the integral representation of the usual SCF 
methods, Hachisu (1986a) took a unique way to impose the 
input parameter for rotation. It uses an axis ratio, r p /r e 
rather than Qc- If the equilibrium stars have no genus, i.e., 
their shapes are spheroidal or quasi-toroidal, r p and r e arc 
the distances to the boundary positions, P along the polar 
axis and Q along the equatorial axis, respectively. Then, 
the ratio lies in the range < r p /r e < 1. For the toroidal 
(donut-shaped) stars, these two distances are the equatorial 
distances to the inner and the outer boundaries respectively, 
and the ratio is designated to have negative values. Actually, 
Hachisu (1986a) found these toroidal-shaped "Dyson- Wong" 
sequences (Dyson, 1983a, b; Wong, 1974) very accurately and 
much more efficiently compared to the previous work by 
Eriguchi & Sugimoto (1981). Note this is not a simple task 
if Qc would have to be chosen as the input parameter, since 
there could be a degeneracy in the rotation speed for a given 
angular momentum (Hachisu, 1986a). 

Our approach is almost the same as the one taken by 
KEH except for having only one gravitational potential to 
solve. The major difference from the Newtonian approach 
is that r e appears explicitly in the equations and should be 
determined during the numerical calculation. If one finds r e 
and Qc, the distribution of the angular speed Q and v can 
be obtained. If the integration constant C is additionally 
known, a final equilibrium solution can be found from equa- 
tion (10) for a given $. These steps are repeated iteratively 
as described below in detail. Note that we recover the uni- 
form rotation if we choose A 2 — > oo while A — corresponds 
to the constant specific angular momentum. 



3.2 Determination of r e ,Q, and C 

To determine r e ,Qc and C, we need to solve equations (9) 
and (10) at three points. Two of them are the two boundary 
positions P and Q. The third one is W where the enthalpy 
or po has its maximum value. Then we have six equations 
for the six unknowns Qp, Qq, Qw, r e , Qc, and C as follows, 

\ In (1 + 2$ P ) + \ In (1 - v%) - \a 2 (Sip - Q c f = C, 



iln(l + 2Sg) + iln(l 



v o) 



\a 2 (q q - n c ) 2 = c, 



lnf/max + i In (1 + 2$> w ) + ^ In (l - Vw) 



- ^A 2 (Q w - Q c ) 2 



A 2 (Q c - ftp) 



C, 

Qprp sin 2 9p 



A 2 (Q C 



Qc 



(1 + 2$p) 2 
Qor\ 



Q 2 r 2 sin 2 Bp 



A 2 (Q c - Qv 



(l + 2<i> Q ) 2 -Q 2 Q r 2 Q ' 
Qwrly 



(l + 2<S> w y 



Q 2 r 2 



(12) 



Since the energy density vanishes at P and Q, we have 
imposed two boundary conditions, H(P) — H(Q) = 1. 
Note that we use p™ ax as a free parameter, so that // max 
is known already. We solve these equations by using the 
Newton- Raphson method. For the spheroidal and the quasi- 
toroidal stars, the number of unknowns is reduced to four 
because two unknowns are predetermined by the relations, 
Qp = Qc and C = \ In (1 + 2$p). Moreover, for the 
uniformly rotating stars the relations can be further re- 
duced, since angular speed is the same at all the positions 
(Qp = Qq = Qw = Qc)- In this case, the values of r e and 
C are easily determined in closed form. 



3.3 Calculation of $ 

In the spherical polar coordinates, the gravitational poten- 
tial $ is given by 



<E>(r) 



-cfr 



s~i I pactivc 

J V^A 

-4ttG / dr I dn 
Jo Jo 

oc 

< ^ /2n(r',r)P 2 „(/i)P2n(/i')pactivc, (13) 



where n = cos 8, P 2 „ is Legendre polynomial of order 2n 
and /2n is defined as 



r /2n + 2 /r 2n+l j f ^ < f 

if r' > r 



fan - \ r 2n J r l2n 



(14) 



It is much more efficient to perform the integration ex- 
pressed with the Legendre polynomials than taking a direct 
3- volume integration. The integration is carried out by using 
the Simpson's formula with second order accuracy. 

Sometimes, we used the overrelaxation technique to im- 
prove the convergence (Varga, 1962). Then, the gravitation 
potential is a weighted average of a newly calculated value 
and the one in the previous iteration step. That is 
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+ (l-wW 



(15) 



where $™ and $ n+ are the gravitational potential at n- 
th and (n + l)-th iteration steps, respectively. The weight 
factor w is determined empirically, so its value is problem- 
dependent. We used the value w = 1 for usual cases (no 
weight) and w — 0.6 for not easily converging problems. 

3.4 Solution procedure 

The computational grid covers from r — to r — y|r e which 
is the same as Hachisu's original one but different from KEH 
who used r max = 2r e . For the initial model of non-rotating 
star before the iteration procedure, we use the modified 
Toleman-Oppenheimer-Volkoff (TOV) equation adapted to 
our configuration. Since this equation is ordinary differen- 
tial equation that depends only on r, we solve it by using 
the fourth order Runge-Kutta method, we used the iterative 
method for the determination of the gravitational potential 
at the center since we do not know the exact value of it. 
After getting the non-rotating model, we decrease the axis 
ratio for the rotating stellar model. The following is a brief 
summery of iteration procedure: 

1) Initial guess of the non-rotating stellar model from mod- 
ified TOV equation 

2) Calculation of the gravitational potential $ from density 
distribution 

3) Determination of the value of equatorial radius (r e ), ro- 
tation velocity at the axis (fie), at the boundary and max- 
imum density position (Op, Oq, flw), and integration con- 
stant (C) using Newton-Raphson method. See equation (12) 

4) Calculation of enthalpy (H) from the value obtained in 
3) and equation (9) 

5) Conversion of enthalpy (H) to the rest mass density (po) 
from equations (5) and (6) 

6) Go back to step 2) until max 

< S, and c "~f G ' J < S, where S is de- 
sired accuracy (iteration criteria). We have chosen 8 to be 
1(T 12 

The calculation of the ring-like toroidal solution starts 
at axis ratio=-0.8. The initial guess of the ring solution 
comes from the Newtonian solution to prevent from the 
failure in getting solutions using Newton-Raphson method 
while solving equation (12), because its solution procedure 
sensitively depends on its initial guess. After calculating ring 
solution with axis ratio=-0.8, we increase the axis ratio to 
get a solution with different axis ratio. 

As for the spheroidal solutions, we start with r p /r e — 1. 
All the procedures are carried out until mass shedding occurs 
when the centrifugal force is so large that the gravitational 
force cannot overcome it. In that case, no stable solution of 
the rotating star exists. 



Table 1. Units of some physical quantities when c = G = M@ 
1 



< 



3.5 Units 

In this paper, we use the units of c = G — Mq = 1. With this 
choice, the length, mass and time units can be automatically 
determined. Consequently, the units of k and po can also be 
determined. Note that the unit of k depends on the choice 



Parameter 


Unit Physical Scale in cgs 


Length 


1.47km 


Time 


4.92 X 10~ 3 ms 


Mass Density 


6.26 x 10 17 g/cm 3 


Frequency (1/Timc) 


2.03 x 10 5 Hz 



of the polytropic index N. In Table 1, we summarize the 
units of these quantities. Our figures are usually expressed 
in these units. 



4 COMPARISON WITH OTHER METHODS 

In this section, we compare our result (pseudo-Newtonian 
Relativistic hydrodynamics approach, pNRH hereafter) with 
three other methods: purely Newtonian, Newtonian rela- 
tivistic hydrodynamics (NRH hereafter) in that special rel- 
ativity is taken into account with Newtonian gravity, and 
general relativistic approaches (GR hereafter). To obtain 
the general relativistic solutions, we use the Whisky _RNSID 
code in Whisky project (http://www.whiskycode.org). The 
code can calculate the equilibrium solutions of rotating stars 
for various equations of states. We have not implemented 
the toroidal configuration in the Whisky code yet. Hence, 
we compare only spheroidal (uniform rotation) and quasi- 
toroidal (differential rotation) stellar solutions. 

For the future reference, we introduce the following 
quantity that measures the importance of the relativistic 
effects, 



R = # max - 1 = (N + 1)k (p^ ax ) 



l/JV 



(16) 



which could be comparable to or larger than unity if the 
relativistic effects are important and R <C 1 for the non- 
relativistic case. Another direct measure of the importance 
of the general relativistic effects is the ration between the 
gravitational potential measured in units of c 2 , i.e., 

GM 

c 2 r eff 

where r e // is the radius of a sphere whose volume is iden- 
tical to the real volume of the star. On the other hand, the 
special relativistic importance due to rapid rotation can be 
measured by the dimensionless parameter 

where J is the angular momentum of the star. Unlike R, 
ip and x should always remain to be less than unity. The 
relativistic effects become significant if these paramegers are 
not very small compared to unity. The angular momentum 
can be computed by 



1 
I 



dSxT^ 
d3xpoH- 



9 

v2 1 



1-d2!1 v a ' 
where g is the determinant of the metric. 



(19) 
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Figure 1. Rest mass density of the uniformly rotating stellar 
models with the axis ratio r p /r e = 0.7 at the equatorial plane 
obtained by four different approaches. Our proposed method of 
pNRH (solid) provides solution that agrees well with the general 
relativistic solution (dashed). The equatorial radius of the New- 
tonian result (dot-dashed) is larger by about 50% compared to 
the general relativistic result as well as our one. For comparison, 
we also have plotted the result from NRH (dotted). Evidently, 
the introduction of the active mass makes further improvement 
of NRH. Note that tp = 0.53 and \ = 0-17 for the model shown 
here, when these quantities arc computed based on our pseudo- 
Newtonian (pNRH) approach. 

4.1 Uniform Rotation 

The density profiles obtained by these four different methods 
for uniformly rotating star with axis ratio of 0.7 along the 
equatorial radius are shown in Figure 1. For this model we 
used N = 1,k= 100, and pS mx = 0.001, so that R = 0.2. 

It is evident that our approach of pNRH gives the den- 
sity profile very similar to that of the GR solution. Our 
pNRH result for the equatorial radius agrees with that of 
general relativity within 5%. On the other hand, the New- 
tonian solution is significantly different from the general rel- 
ativistic one. For example, the the equatorial radius of New- 
tonian solution is about 50 % larger than of the GR method. 

NRH that does not take into account the active mass 
gives the solution closer to the GR one compared to the 
Newtonian approach, but we can see further improvement 
with pNRH. Note the value of R = 0.2 leads to rather small 
relativistic correction. The relative accuracy of our pseudo- 
Newtonian approach over the NRH would be even more ap- 
preciable for models with larger R. Also note that tp = 0.17 
and x = 0.53 for the model shown here, when these quanti- 
ties are computed based on our pseudo-Newtonian (pNRH) 
approach. 

4.2 Differential Rotation 

Figure 2 shows the quasi-toroidal density profiles of differen- 
tially rotating model with r v jr e — 0.35. We took the same 
parameter set with Figure 1 except for P(j lax and A. Since 
the general relativistic code, Whisky _RNSID uses the cen- 
tral rest mass density pg enter as an input parameter rather 
than the p™ ax , we used the value of p c cntcI of the GR solu- 
tion with given p™ 3 ^ as an input parameter for the pNRH. 
The code also uses a rescaled A — A/r e , so we had to find a 



Figure 2. Same as Figure 1 for the differentially rotating star 
with the axis ratio r p /r e = 0.35. The shape of such model is 
quasi-toroidal so that the maximum density does not occur at the 
center (see §4.2 for further details). The similar trend of better 
agreement of pNRH with GR than between Newtonian or NRH 
with GR is evident from this figure. For this model, we found that 
tp = 0.29 and \ = 0.57. 

proper A which could make A 2 = 10 used in the figure. The 
parameters that measure the importance of the relativistic 
effects are found to be tp — 0.29 and \ = 0.57. 

Similar trend with the uniformly rotating model of Fig- 
ure 1 is seen in differentially rotating model of Figure 2: 
progressive improvement from Newtonian approach toward 
the true solution of GR through NRH and pNRH. Our re- 
sult of pNRH (solid) for the density distribution along the 
major axis shows good agreements with the GR solution 
(dashed line). The location of maximum density of pNRH 
differs by 5% compared to that of GR. Newtonian solution 
(dot-dashed) has the equatorial radius which is about two 
times larger than that of GR one. The location of the max- 
imum density is even more different. The difference in these 
quantities between pNRH and GR is much smaller than the 
difference between GR and Newtonian or between GR and 
NRH. Note that the maxumum densities are similar among 
pNRH, NRH and Newtonian, but GR value is slighter larger. 

Not shown in the figures, we have calculated the equilib- 
rium solutions for the Newtonian hydrostatic equation cou- 
pled with the modified Poisson's equation that takes into 
account the active mass. The solutions are not much differ- 
ent from the purely Newtonian case since the active mass 
differs from the rest mass only slightly. We conclude that 
the relativistic treatment of the hydrodynamics is crucial for 
better agreement. Also, it could be improved significantly by 
taking into account the active mass density for relativistic 
cases. 



5 CHARACTERISTICS OF THE 

PSEUDO-NEWTONIAN RELATIVISTIC 
SOLUTIONS 

Newly born neutron stars rotate differentially and are near 
the critical rotation, but turn into uniformly rotating stars 
in relatively short time scale because of the shear viscosity 
and the magnetic tensions. This means that old neutron 
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stars tend to rotate uniformly. We now present the detailed 
characteristics of the pNRH solutions for rotating compact 
stars. We discuss our results separately for uniformly and 
differentially rotating cases. 



5.1 Uniform Rotation 

5.1.1 Spheroidal and Toroidal Solutions 

Examples of our equilibrium solutions for the uniformly ro- 
tating stars are shown in Figure 3. The top and bottom 
panels show a spheroidal and a toroidal solutions, respec- 
tively. We took p™ ax = 0.001 for the spheroidal star in 
the top panel. This choice of p™' 1 * corresponds to physi- 
cal density of a few times larger than the nucleon density 
of 2.7 x 10 14 g/cm 3 . Even though the polytropic equation 
of state is not realistic, these values for k, N, and p™ ax 
are widely used because this set produces a typical size 
and a rest mass density of neutron stars. Our results show 
that the physical size is about 16.6 km which is ~ 5% 
larger than the general relativistic result. With the axis ra- 
tio (r p /r e ) = 0.7, the angular frequency is calculated to be 
597 Hz. For the toroidal star, we take a ten times smaller 
density, Po 1 ^ = 0.0001 which produces a star with a much 
larger size. We note here that the axis ratios for toroidal 
models are denoted by negative numbers following the con- 
vention by Hachisu (1986ab). The axis ratio is selected to be 
-0.35. The size and the rotation frequency obtained are 49.6 
km and 199 Hz, respectively. The critical rotation occurs at 
r p /r e = 0.525 (f ro t = 684Hz) for the spheroidal star and at 
r p /r e = -0.204 (f rot = 228Hz). 

As expected, our solutions show typical features of poly- 
tropic rotating stars. If k increases, the size of star increases 
and the rotation speed decreases. The density falls off more 
rapidly, the size becomes bigger, and the rotation speed be- 
comes smaller as N increases. 



5.1.2 Relation between axis ratio and angular velocity 

In Figures 4, we show equilibrium solution sequences of an- 
gular speed as a function of ellipticity which is defined as 
e = 1 — r p /r e , for different choices of pg^'s, represented by 
different line types, and three differnt values of polytropic 
indices (N), separated by different panels. The values of k 
are differnt for different choices of N: N = 1 and k = 100 
for the top, N = 1.5 and k — 10 for the middle, and N = 3 
and k — 0.1 for the bottom panel. Each panel has three se- 
quences which are three different values of p™ ax : 10~ 2 , 10~ 3 , 
and 10~ 4 . Note that the values of R for the p^ ax = 10~ 2 , 
10~ 3 and 10~ 4 arc 2, 0.2, and 0.02 for N = 1, and R = 1.16, 
0.25 and 0.054 for N = 1.5 and R = 0.086, 0.04 and 0.019 
for N = 3 respectively. 

As shown in the figures, Qc increases with p™ 3 ^ and 
the sequences terminates at the critical rotation speed be- 
yond which the mass shedding occurs. For the Newtonian 
case, due to the simple scaling property of Qc oc VpfJ lax , the 
critical rotation always occurs at the same r p /r e (Hachisu, 
1986a). However, our model shows somewhat different re- 
sults. For the cases of N = 1 and N = 1.5, high density 
models are relativistic but low density models are not. The 
scaling relation of Qc oc v / pQ lax does not hold if the mod- 
els are relativistic. However, we found that the models with 





1e-04 
8e-05 
6e-05 
4e-05 
2e-05 




1 

0.8 
0.6 
0.4 
0.2 




Figure 3. The distribution of po of a spheroidal star (top) and a 
toroidal star (bottom) for k = 100, and N = 1. For the spheroidal 
star, p nax = 0.001 and r p /r e = 0.7 are used as input parameters, 
which are the same for the model shown in Fig. 1. The physical 
size (r e ) and the rotation frequency (f ro t) arc 16.6 km and 597 
Hz, respectively. For the toroidal star, p nax = 0.0001 and r p /r e = 
—0.35. The physical size is 49.6 km and the rotation frequency is 
199 Hz. 



N — 3 are all non-relativistic, and three different models 
with different Qc follows the Newtonian scaling relationship 
very well. The maximum value of ellipticity (or minimum 
values of axis ratio) are also nearly the same for different 

„max 

Po 

Also shown in the inlets of each panel are the 
Qc /Qc,max versus e/emax where the maximum values of Qc 
and e correspond to the critically rotating models. In such 
normalizations, the behavior of Qc with e is almost identical 
for models with different density if the polytropic index is 
fixed. For different N, however, the behaviour of Qc becomes 
somewhat different. Thus we conclude that the scaling re- 
lationship of Qc with e depends on the degree of relativity 
while the functional behaviour is mostly determined by the 
equation of state. 



5.1.3 Constant Rest Mass Sequences 

Most of the neutron stars have their masses near the Chan- 
drasekhar limit (1.4Mq). Even if they undergo dynamical 
changes, the baryon number of the star is conserved. Hence, 
the sequences in fixed rest mass can give an insight into the 
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Figure 4. The rotation angular speed at the rotation axis (fie) 
versus the ellipticity for jV = 1 (T = 2), JV = 1.5 (T = 5/3), 
and A" = 3 (T = 4/3) from top to bottom panels. The three lines 
in each panel represent the results with p™ ax = 0.01, 0.001 and 
0.0001 from the top. The corresponding values of R defined in 
equation (16) from the top are: R = 2, 0.2, and 0.02 for N = 1, 
and 1.16, 0.25, and 0.054 in case oi N = 1.5, and 0.086, 0.040 and 
0.0186, respectively. We also show the normalized angular speed 
versus normalized ellipticity in the inlet of each panel. These fig- 
ures show that the functional behaviour of the angular speed at 
the rotational axis mostly depends on the equation of state. 
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Figure 5. Axis ratio (r p /r e ) - maximum rest mass density (p™ ax ) 
relation in constant masses (1.4Mq and 2Mq) for spheroid (top 
panel) and toroid (bottom panel). 



neutron stars' dynamical properties. The rest mass of the 
star is obtained by 



pQll 



~gd 3 x 



2n dr dOr* sin 9- 



Po 



= (1 + 2$)- 



r , - (20) 

V 1 — v 

In order to obtain the constant mass models, we adjusted 
the maximum density while keeping the equation of states 
unchanged. 

Figure 5 shows the relation between the axis ratio 
(r P /r e ) and maximum rest mass density (po^")- Since more 
elongated figures rotates faster and their sizes are larger, 
they do not need the large value of pcT ax - Accordingly, po lax 
decreases as axis ratio decreases for spheroid and toroid. 
For the toroidal case, the density is much smaller than the 
spheroidal one and its size is very sensitively dependent on 
the value of p^ ax . 

The relation between the rotation frequency at the 
axis versus ratio (r p /r e ) relations is shown in Figure 6 for 
spheroidal (top panel) and toroidal cases (bottom panel). Of 
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Figure 6. Axis ratio {r p /r e ) - rotational frequency relation in 
constant masses (1.4Mq and 2Mq) for spheroid (top panel) and 
toroid (bottom panel). 
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Figure 7. Rotational frequency Q) - T/\W\ relation for the mod- 
els with fixed total masses (1.4Mq and 2Mq) for spheroid (top 
panel) and toroid (bottom panel). 



course more elongated shape and massive star rotates faster 
for spheroidal shape. Rotation frequency of toroid is much 
smaller than the one for spheroid since the specific angular 
momentum is larger for a given rotational frequency. For the 
IAMq spheroidal star, rotational frequency is 539 Hz near 
the critical rotation and 724 Hz for the 2Mq one. For the 
toroidal stars, the critical frequencies become 215 and 280 
Hz for 1.4 and 2.0 Mq, respectively. 

Figure 7 shows the relation between Qc and the ratio 
of the rotational kinetic energy to gravitational potential 
energy (/3 = T/|W|). It is well known that a bar mode in- 
stability occurs if p exceeds certain critical value fid- For 
Newtonian case, the fid for the onset of the bar mode insta- 
bility is 0.2738 (Chandrasekhar , 1969; Shapiro & Teukolsky, 
1983). For the case of fully relativistic stars, the fid lies be- 
tween 0.24 and 0.25 (Shibata, Baumgarte & Shipiro, 2001) 
which is slightly smaller than Newtonian one. The rotational 
kinetic energy (T) and gravitational potential energy (W) 
are defined as 



T = I Jd 3 xnT;^ 

= \j d ^ H Y^2^ (21) 

W = T + M p - M activo , (22) 

where the proper mass M p and active mass M ac ti vc can be 
obtained by 

M P = J d \e^A^ (23 ) 

Mactivc = j d 3 £pactivc\/— 9- (24) 

Like in Figure 6, upper panel and low panel of Figure 7 arc 
for spherpidal and toroidal stars, respectively. We find that 
T/\W\ is much smaller than the critical values of bar mode 
instability for uniformly rotating spheroids with masses 1.4 
and 2.0 M . On the other hand, all the sequences of toroidal 
models have T/\ W\ greater than the instability criteria. Ap- 
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Figure 8. The distribution of po for the quasi-toroidal example. 
The model parameters are the same as that of the sphcrdoidal 
example shown in Figure 3 (i.e., k = 100, N = 1, and p™ ax = 
0.001. The axis ratio of r p /r e = 0.35, and A 2 = 10 are chosen for 
the rotational parameters. The resulting equatorial radius is 17.4 
km, when the baryon mass is assumed to be 1.4 Mq (see Figure 
13). The rotation frequency at the center is fc = 5510Hz but it 
reduces down to 4 % of fc at the equatorial boundary. 



parently, the toriodal configurations are all dynamically un- 
stable. 



5.2 Differentially Rotating Models 

5.2.1 Solution of the Quasi-toroidal Objects 

For differentially rotating stars, the angular frequency at the 
boundary becomes much smaller than the uniformly rotat- 
ing stars in general since the inner parts rotate much faster 
than outer parts. Therefore, solutions could exist well below 
the critical rotational frequency for mass shedding of the 
uniformly rotating stars and the resulting shape of the star 
is quasi-toroidal. In figure 8, we show an example of quasi- 
toroidal distribution of po for the same parameters with the 
uniformly rotating spheroidal star in Figure 3 (i.e., k = 100, 
N = 1, and p™ ax =0.001), except for A 2 which was taken to 
be 10. Also, the axis ratio was chosen to be r p /r e = 0.35. 
The equatorial radius for this model is 17.4 km when the 
mass of the star is assumed to be 1.4 Mq. The resulting 
density distribution is a quasi-toroidal one in the sense that 
density does not vanish at the center. The maximum den- 
sity, however occurs at finite radius so that the location of 
the maximum density is a ring, as already seen from Figure 
2. The rotational frequency at the center is 5510 Hz, but it 
reduces to about 220 Hz at the surface. 

The dependence of the rotation speed distribution on 
the parameter A 2 is shown in Figure 9. For A 2 — 10, the 
rotation speed at the equatorial boundary reduces to 10% of 
the one at the center. As expected, the rotation speed falls 
off more rapidly as A decreases and becomes uniform as A 
increases. 
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Figure 9. The angular frequency distributions of the models with 
the same parameters as that one shown in Figure 8, but three 
different values of A 2 =l, 10 and 100. 



5.2.2 Constant Rest Mass Sequences 

In this section, we describe the characteristics of the models 
with constant baryonic masses of 1.4 and 2 Mq for differen- 
tially rotating objects with ^4 2 =10 and 100. Unlike the uni- 
form rotation, critical rotation for mass shedding does not 
occur in thses values of A. Figure 10 shows that the relation 
between axis ratio (r p /r e ) and maximum rest mass density 
(p™ ax ). Almost similar trend to uniformly rotating models 
is seen but p™ ax falls off more rapidly toward the smaller 
axis ratio. Also the maximum densities are generally larger 
than the toroidal figures of uniformly rotating models, but 
comparable to that of spheroidal models. 

Figure 11 presents relationship between the axis ratio 
(r p /r e ) and the rotational angular speed at the axis of ro- 
tation for the differentially rotating models shown in Figure 
10. Note that the rotational frequency in units of Hz can be 
obtained by multiplying the rotational angular speed shown 
in this figure by 3.23 x 10 4 . Therefore, the rotational freque- 
ncy at the axis of rotation could reach 6130 Hz and 4580 
Hz for 1.4 and 2 Mq models, respectively for the case of 
A 2 — 10. For the case of A 2 = 100, the rotational speed at 
the axis of rotation becomes 1520 and 1170 Hz for 1.4 and 
2 Mq models, respectively. 

Figure 12 shows T/|VF| as a function of the axis-ratio 
for the models shown in Figures 10 and 11. As we stated in 
the previous section, larger A value makes the angular speed 
distribution nearly constant (i.e., uniform rotation). There- 
fore, rotational energy (T) tends to become larger for larger 
A. The most rapidly rotating models with A 2 — 10 have 
T/|VF|=0.17 and 0.21 for the cases with baryon mass 1.4 
and 2 Mq while models with A 2 = 100 can have T/\W\=0.27 
and 0.31 for M = 1.4 and 2 Mq, respectively. Note, for dif- 
ferentially rotating stars, the onset of the dynamical insta- 
bilty depends on the strength of the differential rotation. For 
weak case, it turned out that (3d ~ 0.27 for Newtonian stars, 
just similar to the uniform rotation (see, for example, Liu & 
Lindholm (2001)). However, many stuides have found that 
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Figure 10. Maximum rest mass density (p™ ax ) versus axis ratio 
(r p /r e ) relation for models with two constant masses (1.4Mq and 
2Mq) and two rotational parameters of A 2 = 10 (top panel) and 
A 2 = 100 (bottom panel). 
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Figure 11. Rotation angular speed along rotation axis (He) ver- 
sus axis ratio (r p /r e ) for the models with constant total masses 
(1.4Mq and 2Mq) and differential rotation parameter A 2 = 10 
(top panel) and A 2 = 100 (bottom panel). 



it could be down to fid — 0.2 (Karino & Eriguchi , 2003), 
f3 d = 0.14 (Centrella et al., 2001; Tohline & Hachisu, 1990; 
Pickett, Durisen & Davis, 1996; Liu & Lindholm, 2001; Liu , 
2002), and even down to 0.04 (Shibata, Karino, & Eriguchi, 
2002, 2003). Recent full GR simulation has found (3 d = 0.254 
for weak case (Baiotti et al., 2007). 



6 SUMMARY AND DISCUSSION 

In order to find equilibrium solutions of rapidly rotating 
compact stars which have relativistic motions or relativis- 
tic equations of states, we have proposed an approximate 
method which is appropriate if the gravitational field is not 
very strong. Assuming the weak gravity, the spacetime and 
the hydrostatic equations are derived by only considering 
the Newtonian gravitational potential. However, in order to 
accommodate the relativistic effects, we have adopted the 
active mass density as a source for the gravitational poten- 
tial. The active mass density takes into account all the forms 



of energy density. The numerical calculation has been car- 
ried out by following the Hachisu's SCF method. We have 
obtained the equilibrium solutions for a wide range of pa- 
rameters and topological shapes such as spheroid, toroid, 
and quasi-toroid. Only the polytropic equation of state is 
considered for simplicity. The inclusion of special relativis- 
tic effects could significantly improve the accuracy of the 
solutions if the rotational speed is substantial compared to 
the purely Newtonian hydrodynamical approach. We have 
shown that the solutions can be further improved by taking 
into account the the contribution from the active mass den- 
sity for the computation of the gravitational field. Even for 
a mildly relativistic case (R ~ 0.1), we found that the active 
mass plays important role. 

We have concentrated only on the equilibrium models 
in this paper. The studies on the stability of the relativis- 
tic stars and their physical sequences should be followed by 
integrating dynamical equations. In the subsequent works, 
we will provide the full set of dynamical equations under 
the same assumptions (i.e., Newtonian gravity with active 
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Figure 12. The relationship between T/|VF| and axis ratio for 
stars with M = 1.4M S (solid) and M = 2M Q (dotted) for ro- 
tation parameters of A 2 = 10 (top panel) and A 2 = 100 (bottom 
panel). 



mass and special relativistic hydrodynamics) and carry out 
the dynamical evolution studies of rapidly rotating compact 
stars. 
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